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CHARACTER IZATION  AND  ESTIMATION* 
OF  TWO  DIMENSIONAL  ARMA  MODELS 


R.  L.  Kashyapt 


ABSTRACT 

A  class  of  finite  order  two  dimensional  autoregressive  moving  average  (ARMA)  is 
introduced  having  the  ability  to  represent  any  process  with  rational  spectral  density. 
In  this  model,  the  driving  noise  is  correlated  and  need  not  be  Gaussian.  Currently 

known  classes  of  ARMA  models  or  AR  models  are  shown  to  be  subsets  of  the  above 

”7”A  I  5 

class.  -We  discuss^he  three  definitions  of  markov  property  and  precisely  state^he  class 
of  ARMA  models  having  the  noncausal  and  semicausal  markov  property  without 
imposing  any  specific  boundary  conditions.  Next^ye  consider^the  estimation  of  parame¬ 
ters  of  a  model  to  fit  a  given  image.  Two  approaches  are  considered.  The  first  method 
uses  only  the  empirical  correlations  and  involves  the  solution  of  linear  equations.  The 

second  method  is  the  likelihood  approach.  Since  the  exact  likelihood  function  is 

tu  — 

difficult  to  compute,  resorts^o  approximations  suggested  by  the  toroidal  models. 

The  quality  of  the  two  estimation  schemes  are  compared  via  numerical  experiments. 

K'- 

Finally,  yri  consider^he  problem  of  synthesizing  a  texture  obeying  an  ARMA  model. 


Keywords: 
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I.  Introduction 


Parametric  representations  for  two  dimensional  random  fields  are  useful  in  many 
applications  like  image  synthesis,  classification,  spectral  estimation,  etc.  The  aim  of  the 
paper  is  to  develop  a  finite  stochastic  difference  equation  model  for  regular  two  dimen¬ 
sional  random  fields  having  rational  spectral  densities  and  discuss  related  topics  like  the 
various  definitions  of  weak  markov  processes,  parameter  estimation  and  synthesis  of 
textures  resembling  a  given  non  Gaussian  texture. 

We  will  first  give  the  background  information  regarding  the  structural  representa¬ 
tions.  Rosanov  [l],  Woods  [2]  and  Besag  [3j  have  shown  that  any  Gaussian  markov 
field  having  an  all  pole  spectral  density  (i.e.,  a  reciprocal  of  a  linear  sinusoidal  function) 
possesses  a  finite  difference  equation  representation,  the  so-called  conditional  autore¬ 
gressive  (CAR)  model  in  which  the  driving  input  noise  is  correlated,  but  does  not  have, 
in  general,  a  moving  average  representation.  The  set  of  models  suggested  earlier  such 
as  the  simultaneous  AR  models  (3,  4],  causal  recursive  models  [5,  6]  is  a  proper  subset 
of  the  set  of  CAR  models. 

As  discussed  later,  the  various  types  of  2-D  autoregressive  moving  average 
(ARMA)  models  discussed  in  [7,  8,  9,  11]  have  the  restriction  that  the  denominator  of 
the  corresponding  spectral  density,  say  A(z),  is  factorable,  i.e., 
A(z)  =  D,(z1,z2)D1(zf1,Z2  *)•  Thus  no  general  finite  difference  equation  model  is  avail¬ 
able  for  representing  a  discrete  random  field  having  a  rational  spectral  density  in  which 
both  the  numerator  and  the  denominator  are  not  factorable.  We  emphasize  the  use  of 
the  word  ‘finite’  since  a  simple  spectral  density  such  as 
(1  +  ^(cosXj  +  cosX2)]/[l  ~  0(cosXj  +  cosX2)],  |  <j>\  ,|  0|  <  0.5,  <f>?6,  cannot  be 
represented  by  any  of  the  ARMA  models  in  [7,  8,  9,  11]  using  a  finite  number  of  param¬ 
eters,  but  can  always  be  represented  by  these  models  using  an  infinite  number  of 
parameters.  But  the  principle  of  parsimony  precludes  the  use  of  a  model  having  a  large 
number  of  parameters  especially  in  tasks  such  as  fitting  of  models  to  the  given  data.  In 


contrast  the  class  of  recursive  finite  ARMA  models  in  one  dimension  can  represent  any 
process  with  a  rational  spectral  density  of  finite  order. 

The  2D  case  differs  significantly  from  the  ID  case  in  relation  to  the  markov  pro¬ 
perty.  There  are  3  types  of  weak  markov  property,  namely,  causal  [10],  semi-causal  [8, 
0]  and  non-causal  [1-3].  In  contrast  with  the  ID  case  where  a  process  obeying  an 
ARMA  model  is  a  projection  of  a  vector  markov  process,  the  general  ARMA  model  in 
the  2D  case  is  neither  markovian  according  to  any  of  the  three  definitions  nor  a  projec¬ 
tion  of  another  markov  process.  However,  a  particular  subset  of  ARMA  models  is 
shown  to  possess  the  semi-causal  markov  property  which  was  introduced  in  [8,  9],  We 
will  clarify  the  precise  structure  of  the  ARMA  models  having  the  requisite  semi-causal 
markov  property  without  imposing  any  special  boundary  conditions. 

The  next  topic  to  be  covered  is  the  estimation  of  parameters  in  a  model  to  fit  a 
given  finite  image.  We  present  two  approaches.  In  the  first  approach  the  parameter 
estimates  are  computed  from  the  empirical  correlations  by  solving  linear  equations. 
There  are  no  iterations  in  contrast  with  the  ID  ARMA  model  parameter  estimation 
problems.  The  second  approach  utilizes  the  likelihood.  The  exact  expression  for  the 
likelihood  of  the  given  observations  in  terms  of  the  parameters  is  very  complicated.  We 
consider  an  approximation  which  is  easy  to  handle.  The  approximation  happens  to  be 
the  exact  likelihood  when  the  observations  obey  a  variant  of  the  ARMA  model,  the  so- 
called  toroidal  ARMA  model.  Finally  we  discuss  a  procedure  for  synthesizing  a  texture 
obeying  a  given  ARMA  model. 

Section  II  deals  with  the  general  ARMA  representation,  the  related  markovian  pro¬ 
perties,  and  the  relation  to  existing  2D  difference  equation  models.  Section  III  contains 
the  parameter  estimation  using  the  estimated  correlations.  Section  IV  deals  with  the 
likelihood  approach  which  includes  the  results  of  numerical  experiments  on  the  quality 
of  estimates.  The  next  section  deals  with  the  problems  of  synthesizing  a  texture  to 
resemble  a  real  texture. 


H.  The  ARMA  Model 


Let  y(s),  s  c  L  be  a  two  dimensional  random  field  L  =  {(j,k):  j,k  are  integers}.  Let 
y(*)  be  stationary  and  have  the  correlation  function  Ry(s)  and  spectral  density 
S(z),  z=(z,,z2). 

oo  oo 

s(z)  =  E  E  Ry(j,k)zjz2k,  (2.1) 

j=-oo  k=-oo 

=  t/B(z)/A(z),  z  =  (exp  27riX1(exp27riX2J,  i  =  \ /— T,  (2.2) 

where  A(z)  =  1  -  £)  0rzr,  B(z)  =  1  +  £  Qr  ~  *-*»  =  ^-r  • 

rtNj  rtN2 

N;  are  finite  subsets  of  (L-(0,0))  and  are  symmetric,  i.e.,  if  (j,k)cNj  then  (-j.-kJeN;, 
i  =  1,2.  A  and  B  are  finite  order  polynomials.  A(zj,z2)  =  A(zf*,z2  *).  Similarly  B. 

A(z)  *  0  and  B(z)  *  0  for  all  z  such  that  |  z;|  =  1,  i  =  1,2.  (2.3) 

A(z)  and  B(z)  have  no  common  zero.  (2.4) 

Conditions  (2.3)  and  (2.4)  assure  that  S(z)  is  finite,  and  positive  for  all  real  X  =  (X,,X2). 
Our  aim  is  to  develop  a  finite  difference  equation  representation  for  y(’)  valid  for  any 
spectral  density  S(*)  obeying  (2.2)-(2.4). 

Theorem  1 :  The  stationary  random  field  y(*)  with  its  spectral  density  in  (2.2)  obeying 
the  conditions  (2.3)-(2.4)  obeys  the  bilateral  autoregressive  moving  average  model  which 
can  be  represented  as  in  (2.5)  or  (2.6). 

y(s)  =  E  Ms+r)  +  v^e(s).  (2.5) 

r<N, 
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A(z)y(s)  =  \fv  e(s).  (2.6) 

In  (2.5)  or  (2.6),  e(s)  is  a  zero  mean  stationary  correlated  sequence  with  the  spectral 
density  in  (2.7). 

Se(z)  =  A(z)B(z)  (2.7) 


Proof:  If  y(’)  obeying  (2.5)  exists,  then  taking  spectral  density  of  both  sides  of  the  equa¬ 
tion  (2.6)  and  using  (2.7)  indicates  that  the  spectral  density  of  y(*)  is  as  in  (2.2).  Thus 
all  we  have  to  show  is  the  existence  of  a  sequence  y(*)  obeying  (2.5).  This  will  be  done 
by  construction. 

Let  (w(s),  s  e  L}  be  an  infinite  sequence  of  independent  and  identically  distributed 
random  variables  having  mean  zero  and  variance  unity.  Let  w(*),  e(’),  y(*)  stand 
respectively  for  the  fourier  transform  of  the  infinite  sequence  w(*),  e(*),  and  y(*).  Com¬ 
pute  e(z) 

e(z)  =  \/B(z)A(z)  w(z) 

Then  the  fourier  inverse  of  e(z)  yields  the  sequence  (e(s),  s  c  L}  having  zero  mean  and 
spectral  density  B(z)A(z).  Compute  y(z)  as  shown  below,  which  is  finite  for  all 
|  Z||  =  I  and  |  z2|  =1  in  view  of  (2.3). 
y(z)  =  \fv  e(z)/A(z) 

Rearranging  the  above  equation,  we  get 
(1  “  £  ^r)y(z)  =  v^e(z). 

The  fourier  inverse  of  y(*)  yields  the  desired  sequence  y(*)  obeying  (2.5). 

Comment  1:  In  view  of  (2.7)  the  sequence  e(s)  has  nonzero  auto  correlation  only 
over  a  finite  number  of  lags,  as  displayed  in  (2.8) 
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E(e(t)e(t+s>]  =  -  S  <t>i  *  ^s-r if  s  e  N" , 

=  0,  otherwise 


(2.8) 


where  N;  =  NjU{0,0},  N"  =  {r+s:  r  e  Nj ,  s  e  N2}, 

00,0  =  ^0,0  =  _1) 

<t>8  =  0,  ifsgNj,  9t  =  0  if  r  ^  Nj. 

The  sequence  e(s)  has  non  zero  correlation  with  y(s+r)  only  for  a  finite  number  of 
values  of  r.  To  prove  this  statement,  let  us  find  the  cross  spectral  density  Sey(*)  from 
eq.  (2.6). 

S„M  =  (vtyA(z))  Se(z) 

~  •fv  B(z),  from  (2.7),  (2.9) 

Equating  the  coefficients  of  zr  on  either  side,  we  get 

E(e(s)y(s  +  r)J  =  \fu  <f)T,  if  r  e  Nj,  (2.10) 

=  0,  otherwise. 


Comment  2:  An  alternative  representation  for  y(*)  obeying  (2.5)  is  given  below: 

v/S(z)  y(s)  =  v/S(z)  w(s),  (2.11) 

where  {w(s)}  is  an  independent  and  identically  distributed  (I.I.D.)  sequence  with  zero 


mean  and  unit  variance.  y/A(z)  and  y/B(z)  are  infinite  order  symmetric  polynomials 
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involving  only  a  finite  number  of  parameters  9r  and  <f>r  One  can  verify  that  the  spec¬ 
tral  density  of  y(*)  obeying  (2.11)  is  as  in  (2.2).  The  representation  in  (2.11)  is  more 
suitable  for  synthesis  of  an  image,  as  discussed  in  section  V.  Note  that  the  probability 
density  of  w(*)  can  be  chosen  as  desired. 


Comment  S:  Viewing  (2.6)  as  an  input-output  system  with  e(*)  as  input  and  y(*)  as  the 
output,  we  can  see,  as  in  the  proof  of  theorem  1,  that  a  necessary  and  sufficient  condi¬ 
tion  for  the  BIBO  stability  is  that  A(z)  ^  0  for  |  z,|  =1  and  |  z2|  =1.  In  section  V  a 
specific  algorithm  is  given  for  synthesis  using  this  condition.  The  condition  B(z)^0  for 
|zi|  =1  and  |  z2|  =1  is  needed  for  constructing  a  whitened  representation  of  y(*)  as 
shown  below,  where  w(z)  is  the  fourier  transform  of  the  whitened  sequence  and  simi¬ 
larly  y(-). 


w(z)  =  v/A(z)/B(z)  y(z). 


The  condition  (2.4)  on  A  and  B  in  addition  to  (2.3)  is  needed  to  ensure  the 
identifiability  of  the  parameters  0r  and  <f>r,  as  indicated  later. 


Relation  to  currently  known  ARMA  models 


Case  (1):  The  conditional  auto  regressive  (CAR)  model  [l,  2,  3]  is  a  special  case  of  (2.6) 
and  (2.7)  with  B(z)=l.  The  CAR  models  are  called  as  minimum  variance  representa¬ 
tions  (MVR)  in  [9]. 


Case  (2):  The  simultaneous  AR  model  [3,  9,  10],  also  called  as  a  white  noise  driven 
representation  (WNDR)  in  [9]  is  a  special  case  of  (2.6)  and  (2.7)  where  B(z)=l  and  A(z) 
has  a  factorization  as  in  (2.12). 


(2.12) 


A(z)  —  KDlz)D(z  '1),  D(z)  —  1— 0rzr,  N  need  not  be  symmetric. 

N 

A  simultaneous  ARMA  model  (9,  11,  24]  is  a  special  case  of  (2.6)  and  (2.7)  in  which 
both  A  and  B  have  a  factorization  as  in  (2.12). 


Case  (8):  Consider  the  2-D  ARMA  models  introduced  in  [7]  in  which  4>(z),  a  special  2- 
D  transform  of  the  correlation  function  defined  in  (2.13)  is  a  rational  function  as  in 
(2.14). 


oo  oo 

$(Zl.Z2)  =  S  E  R0lJ2)  zlA> 

i=0  j=-oo 


(2.13) 


=  C(z)/D(z), 


(2.14) 


where 

M,  M2 

C(z)  =  E  E  cij  Azl 

i=0  j=0 


Mj  Mj 

D(z)  =  E  E  du  A4 

i  =0  j  =0 


We  emphasize  that  4>(*)  is  distinct  from  the  spectral  density  S  defined  in  (2.1]  even 
though  $  is  also  called  a  spectral  density  in  [7].  The  ARMA  models  which  possess  a  <!> 
function  as  in  (2.14)  is  a  proper  subset  of  the  processes  having  spectral  density  S  as  in 
(2.2)  and  hence  a  proper  subset  of  the  ARMA  models  defined  in  (2.6)  and  (2.7).  This 
result  is  stated  as  theorem  2. 


Theorem  2:  If  there  exists  a  stationary  process  y(*)  with  its  $  function  as  in  (2.14),  its 
spectral  density  S  has  the  structure  as  in  (2.15). 


S  =  I/B(z)/D,(z)  D,(z— *),  D,(z)  =  E  E  dij  Z1  4 

i=0  j=0 


M2  >  M2,  B  need  not  be  factorable. 


Proof: 


E  E  R(U)  zi  4 


I  — — oo  j=~oo 


oo  oo  .  . 

E  E  R(-i,j)zr  zi 

i=0  j=-oo 


°°  °° 

=  E  E  1  Z2J>  by  replacing  j  by  -j 

i=0  j=-oo 


00  00 

E  E  z2. since  R(*J)  =  RHrj) 

i=0  j  — “oo 


=  $(z,  !,z2  *),  in  view  of  (2.13). 


00 


v  rfn  m\  i  -  G^z2)Gi(z21) 

E  R(0,j)z^  - — , 

j=-oo  G2(z2)G2(z2-«) 


Mji 


where  G;(z2)  =  E  Kikz2 »  *n  v*ew  the  factorability  of  1-D  polynomials 
k=0 


From  (2.1), 


S(zi,z2)  - 


00  00 


0  00 

E  E  +  E  E 

i  =0  j  =-00  i  =-00  j  =-00 


00 


R(iJ)zizi  “  E  R(0>j)z4 

j=-oo 


=  <f(z„z2)  +  ^(zf'jZ^1)  -  Gi(z2)G,(z2,)/G2(z2)G2(z21) 
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_  +  C(z~'l  _  j  G,(z2)|  2 

DM  D(z-')  |  G2(z2)|  2 

_  _ ij3£z] _ 

(D(z)G2(22)1|D(z‘)G2(z^*)| 

=  i/B(z)/D,{z)D,(z~1), 

where  B(z)  is  the  numerator  normalized  so  that  its  constant  term  is  one  and 
D,(z)  =  D(z)G2(z2) 


M,  M2 

Mj2 

x  x  dijzizi 

s  g2jzi 

i=0  j=0 

j=0 

M,  Mj 

=  E  E  dijzizi-  m2  >  m2. 

i=0  j=0 

Comment:  A  simple  consequence  of  theorem  2  is  that  a  process  y(*)  with  spectral  den¬ 
sity  as  in  (2.2)  with  a  non  factorable  denominator  A  cannot  have  a  4>  function  as  in 
(2.15)  and  thus  cannot  have  the  corresponding  ARMA  representation  given  in  [7]. 

Case  (4):  Consider  the  spectral  density  in  (2.2)  in  which  both  A  and  B  have  the  follow¬ 
ing  factorization. 

A(z)  =  KAfzP.fz-1),  B(z)  =  K2D2(z)D2(z_1), 

Dj(z)  =  1  -  £  drzr  ,  D2(z)  =  1  -  £  dr'zr. 
r  t  N4  r  (  Ns 


N4iN5  CL  ,  (0,0)cf  N4t  (0,0)cfN5.  Both  N4  and  N5  are  subsets  of  the  nonsymmetrical 
half  plane  (NSIfP)  L~  indicated  in  figure  1.  Then  the  corresponding  ARMA  difference 


equation  can  be  written  as: 

y(s)  =  E  dry(s  +  r)  +  'fv7  (w(s)  +  E  drw(s  +  r)), 
r<N4  r<Ns 

where  w(*)  is  I.I.D.  (0,1)  sequence.  The  above  equation  is  the  analog  of  the  traditional 
ARMA  model  in  time  series.  If  in  the  above  equation,  in  addition,  B(z)  =  1  or 
equivalently,  d’  =  0,  then  the  corresponding  process  y(‘)  is  said  to  possess  the  weak 
linear  causal  markov  property,  i.e. 

E(y(s) |  all  y(s+r),  r  eL"j  =  £  dry(sfr),  N4  C  L",  (2.16) 

r<N4 

The  corresponding  difference  equation  is  called  as  a  causal  AR  representation  [9]. 
In  this  case,  it  is  possible  to  divide  the  image  at  any  point  s  into  3  parts,  namely,  s  is 
the  present,  the  set  {s+r:  r  e  L-}  is  the  past,  and  {s+r,  r  *  (0,0),  r  £  L~}  is  the 
future. 

< 

Case  5:  Semicausal  models  (9]  is  shown  to  be  a  subclass  of  the  general  ARMA  class  in 
theorem  5  to  be  proved  later. 

Let  us  evaluate  the  conditional  expectation  of  y(s)  given  all  other  values  for  the 
general  model  in  (2.5). 

Theorem  3:  The  sequence  y(*)  in  Th.  1  obeying  (2.5)  and  having  a  Gaussian  density 
has  the  following  conditional  expectation  and  variance 

y,(s)  ^  E[y(s)|  ally(s+r),  r  *  (0,0)}  =  E  Sry(s+r)>  (217) 

r*(0,0) 

EMsJ-y.fs))2]  =*//K,  (2.18) 


where  K  and  gr  are  defined  as 


(2.19) 


K(1  -  E  SS)  =  A(z)/B(z), 
uV 

L'  =  L  -  (0,0),  gr  =  g_r. 

Proof:  Let  G(z)  =  E  gr2'- 
rfL 

u(s)  ^  y(s)  y i(s)  =  (1  -  G(z))y(s).  (2.20) 

The  cross  spectral  density  between  v  and  y  is 
Svy(z)  =  (1  ~  G(z))Syy(l), 

=  i//K,  from  (2.2). 

Hence 

E[u(s)y(s+r)J  =0  Vr  *(0,0). 

Hence  (2.17)  is  true  since  y(*)  is  Gaussian.  To  prove  (2.18), 

Svv(z)  =  (1  -  G(z))2Syy(z),  from  (2.20) 

=  (i//K)(l  -  grzr),  from  (2.2)  and  (2.19). 

Hence  E[v2(s)]  =  i//K. 

XXX 

The  conditional  expectation  in  (2.17)  has,  in  general,  an  infinite  number  of  terms.  The 
next  question  is  the  determination  of  conditions  under  which  the  conditional  expecta¬ 
tion  in  (2.17)  has  a  finite  number  of  terms.  The  answer  is  in  Theorem  4. 


Definition:  A  sequence  y(*)  is  weak  noncausal  markov  if  the  following  is  true: 


E{y(s)|  all  y(s  +  r),  r^0,0]  -  E[y(s)|  all  y(s+r),  rcN,  N  is  finite,  symmetric,  (0,0)  £  N) 


Theorem  4‘  A  stationary  sequence  is  weak  noncausal  markov  and  possesses  a  finite 
linear  conditional  expectation  indicated  in  (2.21)  if  and  only  if  the  process  y(*)  has  an 
all  pole  spectral  density,  i.e.,  B(z)  in  (2.2)  is  a  constant 

EJy(s) |  all  y(s+r),  r^(0,0)]  =  £  gfy(s  +  r),  (2.21) 

«N, 

where  (0,0)  £  N1(  N,  is  symmetric  and  finite,  gr  =  g_r. 

Proof:  ‘ If  ’  part:  Let  the  spectral  density  of  y  be  v/ A(z).  By  theorem  3,  the  conditional 
expectation  is  defined  in  terms  of  gr  in  (2.19), 

gr  =  0r,  if  r  t  Nj, 

=  0,  (0,0),  r  £  N(. 

Hence  the  conditional  expectation  has  a  finite  number  of  terms. 

‘only  if’  part: 

Let  u(s)  =  y(s)  -  y,(s) 

=  y(s)  -  (  E  grzf)y(s)- 

reN, 

Since  y j( •)  is  the  conditional  expectation, 

E[u(s)y(s  +  t)j  =  0,  VI  *  (0,0). 


(2.22) 


(2.23) 


Let  EJu(s)y(s)]  =  K,.  Multiply  (2.22)  on  both  sides  by  y(s  +  t)  and  take  expectation. 
Let  R(t)  =  E[y(s)y(s  +  t)]. 


R(t)  -  E  KrR(t-r) 
r<N 


=  0,  if  t  *  (0,0), 


=  K,,  if  t  =  (0,0),by(2.23). 


Take  fourier  transform  on  both  sides  of  the  above  equation. 


(!  ~  E  Srz)Sy(z)  =  Ki 

reN, 


i.e.,  Sy(z)  has  an  all  pole  spectral  density. 


Comment  1:  Parts  of  the  theorem  4  have  been  known  in  the  literature  [1-3].  The  aim 


of  giving  the  theorem  is  to  show  the  equivalence  of  the  following  three  statements. 


(i)  y(*)  has  the  conditional  expectation  in  (2.21). 

(ii)  y(*)  has  an  all  pole  spectral  density  i//A(z) 

(iii)  y(*)  obeys  the  conditional  AR  model  in  (2.5)  where  the  driving  input  e(*)  has  the 
spectral  density  A(z). 


This  equivalence  is  never  explicitly  stated  in  the  literature.  For  instance  in  [2]  both  (i) 
and  (iii)  are  together  used  in  defining  the  CAR  model. 


Comment  2:  Every  sequence  y(*)  which  is  causal  markov  and  has  the  linear  expectation 
in  (2.16)  defined  by  a  neighbor  set  N4  also  possesses  the  noncausal  markov  property  in 
(2.21)  with  neighbor  set  Nj,  Nj  D  N4.  The  reverse  is  not  true  [3]. 


We  mentioned  earlier  that  y(*)  obeying  a  general  ARMA  model  in  (2.5)  does  not 


possess  the  noncausal  markov  property.  However  a  small  subset  of  ARMA  models 


possesses  another  markov  property  called  as  semi-causal  or  half-plane  markov. 


V  \v V.V.V  W\*“V 


V.\"  V 
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Definition:  (semi-causal  or  half-plane  markov):  y(*)  is  said  to  be  linear  half-plane  mar- 
kov  with  respect  to  the  neighbor  set  N  if 

E[y(s)/  all  y(s+r),  r  e  L+j  =  £  0ry(s+r).  (2.24) 

r<N 

where  L+  =  {(j,k):  k  <  0,  (j,k)  *  (0,0)}.  L+  is  displayed  in  figure  1.  N  is  any  subset 
of  L  +  defined  below. 

N  =  N,  U  N2, 

Nj  =  {(i,0),  (-i,0);  i  =  (2.25) 

N2  {0»j)>  j  — !).•.)  m2,  *  i2,...,±m} 

We  will  presently  display  a  subset  of  ARMA  models,  the  so-called  semi-causal 
models  which  have  the  semi-causal  markov  property. 

A(z)y(s)  =  sTv  e(s)  (2.26) 

where 

A(z)  =  1  -  Yi  Or1*’  N  in  (2.25) 
r<N 

=  i  -  E  +  *rk) 

k  =  l 

m  j  m2 

“EE  ^k,fzikz2f» 
k=-m1  i  =1 

The  correlated  sequence  e(*)  has  zero  mean,  Gaussian  probability  density  and  spectral 
density  in  (2.28)  or  correlation  function  in  (2.29) 


,v  ,*  v\  v,  .%  v%  v.  «*-„  • 

,y.y 


(2.28) 


S,M  *  B(z)  $  1  -  £  *p,o<zf+2fP). 

p=0 

Ree(M)=o  ,  iff*o 

f^ee(^»®)  —  — ^k,0  >  if  k  —  ±l,...,±m„  —  0-j^o  (2.29) 

=  1  ,  if  k  =  0 

=  0  ,  otherwise 

Eq.  (2.26)  can  be  written  as  the  difference  equation  in  (2.30) 

y(i-j)  =  E  My(i+kJ)  +  y(i-k>i)l 

k=i 

+  E  E  0k,t  y(i+k*H)  +  V^e(i,j).  (2.30) 

k=-m, ?  =1 

A  necessary  and  sufficient  condition  for  the  stability  of  (2.30)  is  given  below  [Thm. 
5  of  [9]]. 

A(z,,z2)  /  0  for  |  Zj|  =  1,  |  z2|  >  1.  (2  31) 

The  model  in  (2.26)  is  called  semi-causal  because  it  is  causal  in  the  index  j,  i.e.,  in 
the  RHS  of  (2.30),  j  +  k,  k  >  1  does  not  appear. 

Theorem  5: 


The  stationary  sequence  y(*)  defined  in  (2.26)  and  (2.27)  possesses  the  weak  half¬ 
plane  markov  property  in  (2.24)  if  and  only  if  the  input  sequence  e(*)  in  it  has  the 


correlation  function  in  (2.29)  or  equivalently  y(*)  has  the  spectral  density 
z/B(z)/|  |  A(z)|  |2. 


Proof:  ‘ If  ’  part : 


Sey(z)  =  Cross  spectral  density  of  e(*)  and  y(*), 


^See(z) 

- —  =  i/B(z)/A(z-1). 

Afz"1) 


Expanding  Sey(z)  in  power  series  we  see  that  the  coefficient  of  any  term  involving  z^1, 
0  >  1  is  zero.  Hence, 


R„(M)  A  E|e(i,j)y(i-k,j-f)]  =  0,  iff  >  0. 


R(k)  A  R„(k,0). 


(2.32) 


(2.33) 


Multiply  (2.30)  by  e(i+k,j),  take  expectation  on  both  sides  and  use  (2.32)  and  (2.33). 


R(k)~  £  WR(k  +  P)-R(k"P)l  =ReeM), 

p=l 


(2.34) 


Let  S(zj)  be  the  one  dimensional  discrete  fourier  transform  of  R(k).  Multiply  (2.34)  by 
z~k,  sum  from  k  =  -oo  to  oo  and  use  Ree(k,0)  in  (2.29). 

B(zi)S(z1)  =  i/B(Zl) 

Hence  S(zj)  =  v 


or  Rey(k,0)  =  0,  if  k  *  0 


(2.35) 


Hence  (2.32)  and  (2.35)  yield 
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E(e(s)y(s  +  r)]  =  0  \/r£L  +  ,  (2.36a) 

Since  e(*)  and  y(*)  are  Gauss,  the  above  equation  yields 

E[e(s)|  all  y(s+r),  r  £  L+]  =  0.  (2.36b) 

Taking  conditional  expectation  of  y(s)  given  all  y(s  +  r),  r  £  L+  on  both  sides  of  (2.26) 
and  using  (2.36b)  yields  (2.24). 

‘Only  if’  part: 

Since  (2.24)  is  true,  (2.36)  is  true.  Multiplying  (2.30)  by  e(k,0),  Q  *  j,  and  taking 
expectations  on  both  sides  by  using  (2.36)  we  get 

E[e(i,j)e(k,(? )]  =  0  \jj  *  0  (2.37) 

Multiplying  (2.30)  by  e(i,j)  and  taking  expectation  on  both  sides  using  (2.36)  we  get 

E[e(i,j)y(i,j)}  =  E[e2(i,j)]  =  v  (2.38) 

Multiplying  (2.30)  by  e(i  +  k,j)  and  take  expectations  on  both  sides 

Rey(M)  =  £  Vo(Rey(k  +  P>°)  +  Rey(k-p,0)]  +  Ree(k,0)  (2.39) 

p=l 

By  (2.36),  Rey(k,0)  =  0  if  k  #  0  (2.40) 


Substituting  for  Rey(k,0)  in  (2.39)  from  (2.40)  and  (2.38),  we  get  the  desired  expression 
for  Ree(k,0)  in  (2.29). 
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Comment  1:  For  images  with  specific  boundary  conditions  Jain  [8,  9]  has  shown  that 
the  models  in  (2.26)  have  the  weak  semi-causal  property.  In  theorem  5,  we  have  esta¬ 
blished  the  converse  also  without  imposing  any  specific  boundary  conditions. 

Comment  2:  A  semi-causal  markov  sequence  is  also  causal  markov  only  in  the  degen¬ 
erate  case  B(z)  =  I,  i.e.,  o  =  0  if  k  ^  0.  In  this  degenerate  case,  it  also  possesses  the 
noncausal  markov  property  w.r.t.  to  a  suitable  symmetric  set  N4.  Apart  from  this  case, 
a  semi-causal  markov  sequence  is  never  noncausal  markov  or  vice  versa. 

Comment  3:  In  this  entire  section,  we  have  discussed  only  the  weak  markov  property. 
In  the  Gaussian  case  the  weak  markov  property  is  the  same  as  the  strong  markov  pro¬ 
perty  involving  the  factorization  of  probability  density.  Some  additional  results  con¬ 
nected  with  the  strong  (noncausal)  markov  property  can  be  found  in  [3,  13], 


Ed.  Parameter  Estimation  from  Correlations 


a 
V  \ 


it. 


Given  a  finite  image  over  the  MxM  grid  fl,  we  want  to  develop  a  procedure  for 
fitting  a  ARMA  model  in  (2.5)  to  it,  i.e.,  estimating  the  unknown  parameters  in  it  after 
fixing  the  neighbor  sets  Nj  and  N2.  We  will  give  2  procedures.  In  this  section  the 
parameters  are  estimated  using  the  estimated  correlations  and  it  is  independent  of  the 
density  of  the  image.  The  method  is  computationally  easy,  and  does  not  involve  any 
iteration.  This  state  of  affairs  is  in  contrast  to  the  parameter  estimation  of  ARMA 
models  in  the  time  series  case.  We  will  point  out  the  reason  for  the  difference.  Note 
that  here  0T  =  0(r)  and  Rr  =  R(r).  Nj,  i =1,2,...  are  finite  symmetric  subsets  of 
(L-(0,0)}.  NSi  and  Ng;  are  mutually  exclusive  antisymmetric  subsets  defined  by  the  fol¬ 
lowing  relations: 


..  .V-  -.VV v-.  ..v. v 


N-,  =  NSi  U  N* 


NSi  n  NSi  =  0,  If  (i,j)  6  NSk,(-i,-j)  Cf  Nsk. 


Without  loss  of  generality  we  can  make  NSi  a  subset  of  (L-L~-(0,0))  where  L~  is 
defined  in  Figure  I.  Let  the  polynomials  A  and  B  in  the  spectral  density  of  the  process 
y(')  in  (2.2)  be: 

A(z)  =  l-  £  *r(z'  +  z-')  ,#NSI  =  m,. 
r<Ns, 


B(z)  =  1  +  £  ^r(2'  +  2  ')  #NS2  =  m2- 

r<Ns, 

Let  NSj  —  ~  {S1 . smj) 

The  corresponding  ARMA  model  equation  is 

y(s)  =  S  ^rj)(y(s+ri)+y(s-rj))  +  v^e(s),  (3.1) 

j=l 

where  e(*)  has  the  following  cross  spectral  density 

sey(z)  =  ^  B(z)> 

i.e.,E[e(s)y(s  +  r)]  =  \/u,  if  r  =  0,  (3.2a) 

=  \fv  <f>T,  if  rfN2,  (3.2b) 

=  0,  otherwise.  (3.2c) 

Choose  a  symmetric  set  N3  having  2mj  nearest  neighbors  of  (0,0)  so  that 
N3  D  N2  =  0.  Note  N3  is  not  unique. 

NS3  ~  {tl.t2.  ”.tmi,tfnj  +  1,...,trn)}. 

We  will  obtain  an  explicit  expression  for  the  coefficients  0t ,  <f>t,  v  in  terms  of  the 
correlations  R(s). 


Multiply  (3.1)  by  y(0)  on  both  sides  and  take  expectation  on  both  sides  using  (3.2a). 

R(0)  =  2  £  0(rj)R(rj)  +  v.  (3.3) 

j=l 

Multiply  (3.1)  by  y(s-fsj),  SicNS2,  on  both  sides  and  take  expectation  using  (3.2b). 

R(si)  =  E  0(rj)lR(ri+sj)  +  R(srrj)l  +  ^(si)»  i=l,...,m2.  (3.4) 

j=l 

Multiply  (3.1)  by  y(s  +  t;),  t;£NS3  on  both  sides  and  take  expectation  u:»ng  (3.2c). 

R(*i)  =  E  0(rj)  {R(ti+ri)+R(trrj)>,  i=l,...,m„  (3.5) 

j=l 

In  (3.3),  (3.4)  and  (3.5),  the  true  correlations  R(  )  can  be  replaced  by  their  estimates 
and  the  resulting  equations  can  be  solved  for  0(*),  <f>{‘)  and  v(')  as  indicated  below.  The 
steps  are: 

(i)  Choose  the  set  NS3. 

(ii)  Estimate  the  various  correlation  needed  in  (3.3),  (3.4)  and  (3.5). 

R(r)=  M3J-£y(s)y(s+r), 

1  s 

where  summation  extends  over  all  valid  s  in  Q  and  Mj  is  the  number  of  admissible 
values  of  s. 

(iii)  Solve  the  m,  linear  equations  in  (3.5)  for  0(r|),  .  .  .  ,  0(rm^)  to  yield  the  correspond¬ 
ing  estimates  9 

(iv)  Solve  the  linear  equation  (3.3)  for  v  after  replacing  0(*)  by  9 ,  yielding  the  estimate 


(v)  Solve  the  linear  equation  (3.4)  for  $(•)  after  replacing  6  by  0  and  v  by  u. 

It  is  important  to  note  that  the  computation  does  not  involve  any  iteration.  One 
can  show  that  the  estimates  are  consistent,  i.e.,  as  the  size  of  the  image  M  goes  to 
infinity,  the  estimates  tend  to  their  true  values  provided  A(z)  and  B(z)  do  not  have  any 
common  factors. 

Comment  1:  The  procedure  needs  the  choice  of  the  set  N3.  The  set  could  be  arbitrary 
as  long  as  it  is  exclusive  of  N2.  Empirical  evidence  indicates  that  the  one  suggested 
here,  namely  having  2m,  nearest  neighbors,  leads  to  the  estimate  with  higher  accuracy 
than  other  choices. 

Comment  2:  As  noted  earlier,  no  iteration  is  needed  in  the  computation.  In  contrast, 
the  estimation  of  parameters  by  the  covariance  method  in  the  one  dimensional  ARMA 
model  is  much  more  complicated  and  involves  iteration.  The  reason  for  the  different 
behavior  is  the  difference  in  the  model  equations.  In  the  2D  case  the  input  sequence  e(’) 
is  correlated.  In  the  1-D  ARMA  case  the  input  sequence,  w(*),  is  independent.  If  we 
convert  the  1-D  ARMA  equation  into  a  form  similar  to  eq.  (3.1),  then  the  computation 
procedure  indicated  in  this  section  can  be  used  for  1-D  case  also.  Note  that  in  this  sec¬ 
tion,  we  aim  at  estimating  directly  the  coefficients  occurring  in  the  spectral  density, 
whereas  in  traditional  1-D  ARMA  case,  we  estimate  the  coefficients  of  C(z,)  and  D(z,) 
where  S(z,)  =  |  |  C(z,)|  |  2/|  |  D(zi)|  | 2. 

The  sequence  of  computation  is  illustrated  by  an  example. 


Example  1:  Let  the  ARMA  model  be  as  in  (3.6)  and  (3.7). 

y(s)  =  e  £  y(s+r)  +  Vu  e(s)  , 

rtN, 


Nj  =  {(0,1),  (0-1),  (1,0),  (-1,0)} 

B(z)  =  1+  ^E*'N2  =  N„ 

r(Nj 

The  required  choice  of  N§3  =  (1,1)-  Note  R;  j  =  R;_j  =  Rj ;.  From  (3.4),  we  get 

Ro,o  =  ^Rli0  +  v- 


From  (3.5),  we  get 


Ri,o  ~  0[Ro,o  +  R2,o  +  2Ri,i]  + 


From  (3.6),  we  get 


Rl,l  -  2^[R1,0  +  Rl,2l- 


(3.10) 


Solving  (3.8)-(3.10)  for  0,  <f>  and  v  is  straightforward.  The  numerical  results  are  given  in 
the  next  section. 


IV.  Likelihood 

When  the  number  of  parameters  to  be  estimated  is  not  very  small  compared  to  the 
image  size  M2,  then  the  estimates  given  earlier  may  not  be  accurate.  Hence,  we  intro¬ 
duce  the  more  accurate  method  of  estimation,  the  so  called  likelihood  method.  As 
before  let  the  given  set  of  observations  be^. 

v  =  Col.[y(s),s«n) 
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n  =  {(*ij))0  <  i,  j  <  M-l} 


Let  the  spectral  density  of  y  be  t/B(z)/A(z).  Let  us  assume  that  y  is  Gauss.  Then  we 
can  find  the  correlation  function  R(s)  of  y(*)  as  a  function  of  0;,^;  and  v  from  the  spec¬ 
tral  density.  Then  we  can  find  the  covariance  of  the  vector  £,  say  C.  Thus  ^  is  Gauss 
(O,C(0,^,t/)].  But  the  matrix  C  is  of  dimension  M2xM2  and  R(s)  is  not  a  simple  function 
of  s,  0,  v  and  <j>.  Hence  the  above  density  expression  for  ^  is  not  useful  for  tasks  like 
maximizing  it  to  find  the  parameter  estimates.  We  have  to  be  content  with  an  approx¬ 
imation  to  the  probability  density  function  of^  so  that  it  is  amenable  for  optimization. 

Let  {Yr,rcQ}  be  the  finite  fourier  transform  of  the  finite  sequence  {y(s),sffi}.  Then 
as  M  tends  to  infinity,  the  sequence  {Yr}  is  independent  and  Gauss  with  mean  zero  and 
variance  m2  M  in  (4-2)- 


SrLW  = 


A(z  =  exp(i(27r/M)r)) 


•,  i  = 


1-0tQj 

r\J 


=  Col.(2cos(27r/M)rTs,  scNS2),  «r  =  Col.(2cos(27r/M)rTs,  stNSl) 


&  =  Co\\<t>v  rcNS2))r0  =  Col.(0r,  reNg,) 

Asymptotically,  the  probability  density  of  {Yr}  is: 

p(Y„  =*  [  n  l/(2*S,(«,£p)l'/2exp[-(l/2)  £||Y,||  */M*S^)].  (4.3) 

ref)  r(f2 

We  can  show  [12]  that,  if  we  transform  {Yr,  rcfl)  in  (4.3)  into  {y(s),  scfl},  then  the 
RIIS  of  (4.3)  is  the  exact  probability  density  of  a  set  of  observations  obeying  the 
toroidal  variant  of  the  ARMA  model  described  below  which  is  valid  for  the  region  fi 
only. 


(4.4) 


■I’'.  'y,’'T’Y 
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y(s)  =  E  ^r(y(s  ©  r)  +  y(S  ©  ^)  +  y/v  e(s),  sen, 

reNg! 

where  ©  denotes  summation  modulo  M, 

y(s)  —  y(s  mod  M),  e(s)  =  e(s  mod  M). 

fe(s),  seH}  have  the  correlation  function  described  earlier. 

The  expression  on  the  RHS  of  (4.3)  has  to  be  maximized  w.r.t.  £  and  v.  It  is 
more  convenient  to  work  with  the  log  likelihood. 

Jit&&1')  =  ~(M2/2)ln(2^)-l/2  £log[l  +  ^TV>r)/l-0Tar)] 

r 

-(1/2)  E 1 1  Y,I  i2(i-eT<.r)/K(i  +4H,)yP,  (4.5) 

ren 

Maximizing  Ji  w.r.t.  v  yields 

V  =  (El  I  Y,|  I  2( !-/»)/( 1  +  A))/M<  (4.6) 

r 

Substituting  it  back,  and  simplifying,  we  see  that  the  ML  estimates  of  £  and  ^  are 
obtained  by  minimizing  J(0,^>)  w.r.t. £  and 

HiS  =  EMU  i-iTa)l 

reO 

Since  the  minimizing  value  of  [0,<f> )  has  to  yield  a  finite  value  for  J(-),  the  ML  estimates 
of  6  and  <f>  automatically  satisfy  the  conditions  and  1  +^>T^r  J*5  0  lHp'o,.  /  0  for  all  r. 
Thus  the  ML  estimates  of^  and  ^  satisfy  (2.3).  We  cannot  make  such  a  claim  for  the 
estimates  obtained  by  the  correlation  methods  especially  for  small  M.  The  numerical 
aspects  of  maximization  have  been  discussed  in  (4,15]  for  the  case  of  AR  models. 


The  likelihood  approach  can  be  adapted  to  the  particular  situation  on  hand.  If  we 
know  that  the  observation  is  the  sum  of  a  signal  obeying  a  CAR  model  and  an  additive 
noise,  then  we  can  directly  write  the  likelihood  and  estimate  the  parameters  of  the 
CAR  model  and  the  corresponding  spectrum.  If  the  signal  plus  noise  assumption  is 
true,  any  spectral  estimation  method  which  ignores  the  noise  will  not  give  good  results. 
This  feature  has  been  documented  in  [16]. 

Numerical  Experiments 

The  correlation  and  maximum  likelihood  estimates  are  compared  via  numerical 
experiments.  The  image  model  in  example  1  of  section  IV  is  considered,  with  numerical 
values  9  —  .22,  <j>  =  .2  and  v  —  1.  Ten  different  images  of  size  64x64  obeying  this 
model  were  synthesized  using  different  random  sequences,  as  discussed  in  section  V.  In 
each  case,  the  parameter  estimates  were  computed  by  both  the  methods.  For  <j>,  9  and 
v ,  the  mean  of  the  10  estimates,  the  standard  duration  (SD)  of  the  10  estimates  and  the 
root  mean  square  value  of  the  deviation  of  the  estimate  from  the  true  value  (RMSE) 
are  given  in  table  1  for  both  correlation  and  ML  estimates.  Similar  experiments  were 
performed  with  images  of  size  32x32  and  16x16  and  the  results  are  also  given  in  table 
1. 

For  64x64  images,  the  SD  and  the  RMSE  are  close  to  one  another.  Further  the 
correlation  estimates  and  the  ML  estimates  of  <j>  and  9  have  similar  accuracy.  The 
correlation  estimate  of  v  appears  to  be  slightly  biased.  But  as  the  size  of  the  image 
decreases,  the  RMSE  values  of  the  ML  estimates  are  less  than  the  corresponding  values 
of  the  correlation  estimate.  This  feature  is  to  be  expected.  But  the  quality  of  correla¬ 
tion  estimate  is  not  unduly  low.  For  instance  for  32x32  image,  the  RMSE  values  for 
ML  and  correlation  estimates  of  9  are  .0124  and  .0182,  not  very  drastic.  In  many  image 
processing  problems  the  correlation  estimates  appear  to  be  adequate,  especially  in  view 
of  their  low  computational  demand. 
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V.  Synthesis 

An  interesting  problem  in  image  processing  is  the  synthesis  of  an  image  which 
resembles  a  real  texture.  There  are  many  methods  of  synthesizing  images,  each  one 
based  on  a  different  type  of  model.  Synthesis  has  been  done  using  various  types  of  2D 
AR  models  [5,  13,  17-20J,  and  mosaic  models  [21-22).  Our  aim  is  to  explore  the  use  of 
ARMA  models.  In  an  ARMA  model  in  (2.11),  there  are  three  sets  of  parameters: 

(i)  the  order  of  ARMA  model,  i.e.,  the  sets  NSj  and  NS2 

(ii)  the  values  of  the  coefficients 

(iii)  the  histogram  or  density  of  the  input  process  w(*)  in  (2.11). 

The  choice  of  the  appropriate  order  of  the  2D  AR  models,  suitable  for  a  given 
image  is  given  in  [12,  23).  A  similar  procedure  is  suitable  for  ARMA  model.  The  esti¬ 
mation  of  parameters  has  already  been  discussed.  We  will  discuss  a  convenient  method 
of  synthesizing  an  image  obeying  a  ARMA  model.  Techniques  for  synthesizing  images 
via  CAR  and  SAR  models  to  resemble  a  given  texture  which  takes  into  account  all  the 
aspects  mentioned  above  is  given  in  [13,  19). 

The  synthesis  of  a  finite  MxM  image  to  obey  exactly  the  difference  equation  in 
(2.5)  or  (2.11)  is  very  difficult.  It  involves  the  factorization  of  a  M2xM2  matrix  whose 
elements  are  the  correlation  of  various  lags.  It  is  still  more  difficult  to  ensure  that  the 
density  of  the  synthesized  y(*)  has  the  prespecified  form.  Instead  consider  the  toroidal 
variant  of  the  ARMA  model  in  (2.11)  which  can  be  compactly  written  in  terms  of 
{Yr,  rcH},  the  FFT  of  (y(s),  scU) 

V^Yr  =  v^Wr,  (6.1) 

Ar  =  A(z  =  exp(\/-T(2ff/M)r),  Br  =  B(z  =  expv/-l27rr/M). 

{Wr,  reV}  is  the  FFT  of  (w(s),  scfl},  (w(s),  seQ}  being  drawn  from  an  I.I.D.  sequence 
with  zero  mean,  unit  variance  and  the  histogram  P  mentioned  later  on.  As  M  tends  to 
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infinity,  the  second  order  properties  of  the  sequence  obeying  the  toroidal  ARMA  model 
in  (6.1)  tend  to  that  of  the  general  ARMA  model  in  (2.5)  or  (2.11).  Eq.  (4.4)  is  also  an 
equivalent  toroidal  representation  which  is  the  direct  analog  of  (2.5). 

To  generate  the  histogram  P  of  w(*),  for  the  given  image  we  proceed  as  follows. 
Using  the  given  image  say  (y'(s),  seQ}  and  the  (estimated)  parameters  6,  <f>,  v  etc.,  gen¬ 
erate  residuals  {w*  (s),  stfi}.  The  corresponding  FFT  can  be  computed  as 

w;  =  y;  -JxJbTv 

> 

where  {Y,,  reft}  is  FFT  of  (y'  (s),  sell}.  The  inverse  FFT  of  {Wr,  rcfl}  yields  the 
(w'(s),  sefl}.  The  histogram  of  v!  is  the  required  histogram  P. 


The  Synthesis  Procedure  is  as  follows: 

(i)  Generate  a  sequence  {w(s),  scfl}  drawn  from  an  I.I.D.  population  with  zero 
mean,  unit  variance  and  histogram  P. 


(ii)  Compute  Wr,  refi 


wr  =  £ 

sell 


where  fr!I  is  the  fourier  array,  r<ft,  seQ 

fr,s  =  expjx/-!  sTr),  s,refl 

(iii)  Compute  ^(s),  scfl  the  required  image  matrix 


.  -  /.y*  •  .  «  .  *  s  •  i.  '  ^  •  .  •  » *  ■  V  V  V  ' 


VA’.'V 

v'Av 


,-;.v 


VI.  Conclusions 


We  have  introduced  the  general  class  of  two  dimensional  ARMA  models  which  can 
represent  any  discrete  rational  spectral  density  and  shown  that  the  various  classes  of 
two  dimensional  difference  equation  models  discussed  in  the  literature  are  subclasses  of 
this  general  class.  We  have  also  given  various  definitions  of  weak  markov  processes  and 
precisely  characterized  subclasses  of  ARMA  models  having  the  various  types  of  markov 
properties.  Two  methods  are  given  for  estimating  the  parameters  in  the  model. 
Finally,  a  technique  is  given  for  synthesizing  an  image  obeying  a  given  ARMA  model. 

Acknowledgement:  The  author  would  like  to  thank  Mr.  G.  Boray  for  carrying  out  the 
numerical  experiments  of  section  IV. 
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